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UNSTEADY TRANSONIC SMALL DISTURBANCE THEORY V7ITH 
STRONG SHOCK WAVES 


1 . INTRODUCTION 


The most common methods of predicting steady flow aero- 
dynamic characteristics at transonic speeds are either the 
Transonic Small Disturbance (T3D) theory (ref. 1) or the Full 
Potential Equation (FPE) theory (ref. 2) . The more accurate 
Euler equations solutions (ref. 3) are expensive to obtain, 
although for flows with strong shock waves such solutions are 
essential. The FPE theory is based on the assumption that the 
flow is isentropic and irrotational and generally has a (numeri- 
cally) exact treatment of the wing boundary conditions. The 
TSD theory is an approximation to the FPE theory and ±hin wing 
boundary conditions are used in the solution procedure. One of 
the advantages of the TSD theory is the flexibility in deriving 
the approximate equation. This flexibility is generally 
utilized by a choice of a transonic scale parameter. The basic 
assumptions of isentropy and irrotationality in both these 
theories are only valid v;hen there are no shock waves in the flov/ 
or when any shock waves are weak. The generally accepted defini- 
tion of a weak shock is v;hen the local Mach number just ahead 
of the shock is less than 1.3. When both TSD and FPE solutions 
are compared to the more realistic Euler equation solutions it 
is found that the agreement is satisfactory provided that the 
basic restriction to weak shock waves is not violated. The use 
of thin wing boundary conditions can also introduce errors into 
the TSD solutions. If the flow has strong shock waves, however, 
then there is considerable disagreement between both TSD and FPE 
solutions and Euler equation solutions. Generally the predicted 
shock location for the potential theories is much further aft 
than that for the Euler equation solutions. This is due to the 
isentropic assumption being invalid in these flows. The causes 
of the error in the shock location in the steady TSD theory for 
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two-dimensional flow have been examined in reference 4 where a 
correction procedure has been derived that allows the basic 
equation within the formal accuracy bounds of a small disturbance 
theory. The basic hypothesis of the theory is that the error 
in shock location is primarily due to the stronger shock 
predicted by TSD theory compared to the shock strength of the 
Euler equations. It is also assumed that if the shock strength 
is suitably corrected then the shock location should be approxi- 
mately correct. The technique used in reference 4 makes use of 
two TSD solutions, with different scaling parameters, and an 
interpolation scheme derived for discontinuous transonic flows. 
Examples of flows with strong shocks computed with this method 
agree satisfactorily with the Euler equation solutions, although 
the use of the thin airfoil boundary conditions in the TSD theory 
can give rise to errors near the leading edge. 

The aim of the present paper is to extend the basic concept 
of reference 4 to include unsteady transonic potential flows. The 
only satisfactory unsteady transonic method available is the low 
frequency theory of Ballhaus and Goorjian (ref. 5) which numeri- 
cally integrates the nonlinear low frequency transonic small 
disturbance equation in a time accurate manner. As in steady flow 
the results for a small disturbance equation formulation compare 
satisfactorily with solutions of the Euler equations when any 
shock waves in the flow are weak. However, the accuracy of the 
solution dxminishes, particularly as regards the shock location 
and motion, as the shocks become stronger. Again, as in the 
steady flow, it is assumed that this error is due to the shock 
strength of the small disturbance solution being larger than the 
corresponding Euler equation solution and that it is this dif- 
ference that leads to the wrong shock location. 

In reference 4 the correction to the TSD equation is obtained 
by computing two steady state solutions and then using an 
interpolation technique to give the required solution. This 
technique is not really feasible for unsteady flow since the 
correction procedure is required for each time step in the 
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solution. As there may be 1,000 or more such time steps the 
technique used in reference 4 would require an inordinate amount 
of computer time to obtain a solution. Consequently, a different 
means of implementing the shock strength correction is developed 
in this paper. The technique involves the addition of higher 
order terms, which are formally of negligible magnitude, to the 
low frequency TSD equation. These terras are then chosen such that 
any shock waves in the flow have strengths approximately equal 
to the appropriate Rankine-Hugoniot shock strength. Two correcting 
approaches are investigated in the paper. The first is to 
derive a correction for the mean steady flow and then simply use 
this corrected form for oscillatory flows. The second is to 
derive a correction for both steady and oscillatory parts of the 
flow. This second development is the most satisfactory and com- 
parisons of the present results with Euler equation results are 
generally favorable, particularly regarding shock location, 
although there are some discrepancies in the pressure distribu- 
tion in the leading edge region. 


2. BASIC SMALL DISTURBANCE FORMULATION 

The low frequency TSD equation for the perturbation velocity 
potential 4>(x,y) at a free-stream Mach number M^ can be written 
in the form (ref. 6) 

5 2M^c 

(1 - M )(J> + (j) = (y + DM'^cI) (J) „ + 6 ^ (1) 

oo' ^xx ^yy ' oo^x XX U ^xt ' ' 

^OO 


where y is the ratio of specific heats and q is the transonic 
scaling parameter. The two most commonly used (ref. 6) values of 
q are 2 (Spreiter scaling) and 1.75 (Krupp scaling) . In the 
transonic limit of M^ -»■ 1 both scalings are identical. The 
pressure coefficient Cp(x,y,t) is given in the low frequency 
small disturbance theory by 



Associated with Equation (1) are the usual tangency and far 
field boundary conditions. The weak shock jump conditions for 
Equation (1) are 

2 2M^cx 

( 3 ) 


where jj^ Jj denotes a jump through a shock wave, (}>g is the angle 
that the normal to the shock makes with the x-axis and 
Xg = shock speed. For a shock normal to the free 

stream the shock strength, a^, is defined as 


a^(t) = Cp(t) - Cp(t) = -2 


Cl(t) - C 


4M^cXg (t) 

P (Y+1)M% J 


(4) 


where C^, C~ are the pressure coefficients just ahead of and 
behind the shock and 


-2(1 - M^) 

C* = (5) 

P (Y + 1 )M^ 

For a harmonic oscillation with frequency w 

X = 1 6 x joje^^^ ( 6 ) 

where | ] is the magnitude of the shock motion. Now 
|5Xg| ~ \/here 63 ^ is the amplitude of the oscillatory 

motion and generally 163 ^! << 1. Since in the low frequency 
equation. Equation (1) , it is implied (ref. 6 ) that 
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k = ^ ~ O(B^) « 1 


it can be seen from Equations (4) and (5) that 


4M^cx 
00 s 


„ - 0(6,k) << C* 

(Y+1)M% ^ P 


Hence, the shock speed can be neglected in the computation of 
the oscillatory shock strength. Thus 


= -2{C+(t) - C*} (7) 

For the following shock jump relations for the Euler 
equations it is assumed that the low frequency assumption noted 
above still applies. That is, the shock speed is sufficiently 
small that the quasi-steady relation holds. 

Consider now the quasi-steady Euler equation normal shock 
jump, ag(t), in terms of and C'*'{t) which is given (ref. 7) 
by ^ 






ym 


o + C^{t) 
2 p 


( 8 ) 


where the upstream shock Mach number is given by 


V 


Mp(t) = f[2 + (Y - DmJ]/ MV(t) + l] ^ ^ - 2 ) 

® (Y - 1) CO ~ p J 


(Y - 1) 


(9) 


The TSD (with q = 2.0) and Euler shock strengths are plotted 
against for M^ = 0.755 in figure 1 and it can be seen that as 
|Cp| increases the discrepancy between o,j, and increases. It 
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should be noted here that different transonic scalings not only 
give a different value of C* but generally a different value of 
Cp. Thus, for different scalings the shock strength may vary 
considerably . 

The basic hypot sis of the present theory, as in reference 
4, is that the error in the shock location in both the steady 
and unsteady TSD solutions is due primarily to the error in the 
shock strength as exhibited in figure 1. Solutions of the Euler 
equations give a different pressure jump across a shock wave 
than potential theories because of the entropy production due 
to the shock wave. It is suggested that if th^ TSD equation is 
altered, still within its formal accuracy bounds, such that the 
shock jump approximates the Euler equation shock jump, then the 
resulting equation is a better compromise in representing the 
actual flow. The reason for this statement is that by matching 
the shock jump the new equation implicitly introduces an addi- 
tional mechanism (formally negligible) that cancels the entropy 
production and rotationality errors in a potential formulation. 

The method introduced in reference 4 to modify the TSD 
equation is to interpolate two TSD results with different scalings, 
to obtain the desired pressure jump. While this technique is 
adequate for steady flow corrections it is not really adequate 
for unsteady flow results since it would require two TSD results 
at each time step. This would lead to an inordinate amount of 
computer time. Consequently simpler means of modifying the TSD 
equation are examined; these are discussed in the next section. 


3. MODIFIED SMALL DISTURBANCE EQUATION 

In an initial attempt to derive a small disturbance equation 
that will give a more realistic shock jump the following equation 
is used. 


(a + b(J)^ + 


2 

Cd ) d) "t (b 

^xx ^yy 




xt 


(10) 
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If a = 1 - M^, b = (y + 1)M^ and c = 0 Equation (10) is identical 
to Equation (1) . 

In the following analysis it is assumed that all shock 
waves are normal to the free stream. In this case contributions 
from the tern to the jump relation can be neglected. 

Consider for the moment the more general form of Equation (10) 

+4) -2M^ — <1) ^ = 0 (11) 

^xx ^yy «> n xt 

00 

The normal shock jump relation is 


F(oJ,o) - 2mI X = 0 

X y S 

**^co 


(12) 


where 4> is the value of 4> just ahead of the shock and a is the 
shock strength. The function F(<)) ,o) is given by 

X 


F(4^,a) = I 






(13) 


. + o 
" 2 


The problem is to pick a suitable form of f (tj) ) , and in 
the following discussion some conditions that should be satisfied 
by suggested. 

(a) In general the function f(c}) ) will be nonlinear in 

•X 

which leads to the possibility of multiple parabolic points 
when 


f(<l>j^) = 0 (14) 


It is advisable that at least one of the roots of Equation (14) 

be (f) which is the parabolic point for a conventional (say Krupp) 
*c 

TSD solution. Thus 
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f(0* ) = 0 


(15) 


In order to avoid unrealistic multiple parabolic points in the 

domain of interest it is also desirable that f(c() ) be a monotoni- 

cally decreasing function in some range 4^ • < < v>v • Thus 

J ^ ^ ^^min — — ^^max 


df 

d(o^) 


< 


0 for - “^x 

min max 


(16) 


(b) The normal shock strength should be the same as the 
normal Euler shock strength. Neglecting the shock speed term 
in Equation (11) , as discussed in the previous section, this 
gives 


F(=;^, cg) = 0 


(17) 


where is given by Equation (8) . 

These three conditions, given in Equations (15) , (16) , and 
(17) , will give the desired shock strength with no undesirable 
side effects. 

Equations (15) , (16) , and (17) must be satisfied at each 

iterative step or time step in order to get the correct shock 
jump. While there is little option to satisfying each of these 
equations at all iterative steps at this preliminary stage of 
the investigation, it is possible to reduce the amount of com- 
putational work required for an unsteady example if a steady 
state solution is first computed. 

It is ofren the case that the unsteady pressure distribution 
is related by a small steady perturbation to a mean state, which 
itself is close to the steady state result. Hence, if a modified 
small disturbance equation has been evolved for the steady state 
case it may be possible to use a simple analytic perturbation 
of the function f(<^„) that will give good accuracy in the 
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neighborhood of the steady state result. The additional condi- 
tion for this idea is given below. 

(c) For the shock strength to be approximately correct 
over a range of values close to the steady state value of 
<{),(}) , Equation (17) can be expanded as a Taylor's series to 

give 


'■'♦x ' "e * 

S S 



’ 3F ' 


^'^eI 


1 

+ 

^x_ 




(18) 


where the subscript s denotes a value in the steady state. 
Now by definition 


F(4^ , Og ) = 0 

s s 


(19) 


and hence for Equation (18) to be satisfied for a range of 
'♦x - ' 


5F 3F 


0 


( 20 ) 


For an unsteady flow therefore there are two options for practical 
calculations . 

Option I.- The function f(^„) can be chosen such that 
Equations (15) -(17) are satisfied at each iterative and time 
step. 

Option II.- Equations (15), (16), and (19) are satisfied 
for each iterative stop in the steady calculations and Equations 
(15) , (16) , and (20) are satisfied at each time step for the 

oscillatory calculation. 
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4. PRELIMINARY IMPLEMENTATION OF THE CORRECTION PROCEDURE 

It is possible to derix'^e functions f (tj) ) that satisfy all 
of the conditions given in the preceding section, but in a limited 
study, such as the present one, it is more convenient to satisfy 
only certain of these conditions explicitly and to test the 
resulting function with respect to the other conditions. In the 
following discussion the set of conditions given in Option II 
are considered. 

The most crucial conditions to satisfy are the shock 
strength conditions. Equations (19) and (20) since the object of 
the present excercise is to realize the correct shock strengths. 
Consequently, these conditions are satisfied explicitly. For 
simplicity in derivation, it is assumed that the modifications 
to the small disturbance equation will be sufficiently small that 
Equations (15) and (16) V7ill be, at least, approximately 
satisfied. This aspect v/ill be considered with each application 
of the theory . 

The form of f(<J>jj) chosen is that of Equation (10), that is 


f((i),^) = a + b(|)^ + C(J)J (21) 

where a, b, c are either constants of simple functions of 

The particular form of Equations (15) , (16) , (19) , and (20) 

is now 


^ * *2 
a + b^.^ + c4>^ = 0 


( 22 ) 


b + 2cd>^ < 0 for ^ < <l>j^ < 

min max 


f cr 1 



5+ - 

1 "a = J 

+ c 



(23) 


a + b 


0 


(24) 
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aE 



1 - ® 

-4- C* 

Jj 

^ 2 




1 . 


s s s s s s 


■ da ■ 


f \ 

. + *"^s 

- , ■V 

+ 

‘f’x - 2 

V J 

s 

s 

V, y 


db ' 

+ 

dd)'^ 

V 

s 


2 1 


+ + ^■ 
6^ - a„ d) + — :5-^ 

E ^x 3 

s s s 


dc 


\ J 


= 0 


(25) 


Two applications of this correction were developed as 
follows. 

Method I.- The coefficients b and c were chosen such that 

2 

for the mean steady case Equation (24) is satisfied with a = 1 - 

and b = (y + 1)M*^. The coefficients are then frozen for the 

00 

unsteady part of the computation. Thus, no further correction 
is required for the unsteady calculation. 

Method II.- The coefficient c is chosen such that Equation 
(24) is satisfed and 

a = 1 - M^ 

oo 

b = (y + 

The exponent a is chosen such that equation (25) is satisifed. 

This form of f (4>v) does not have the same propensity for 
producing multiple parabolic points, as altering b for the steady 
state solution but incurs the cost of altering f (0^^) at each 
time step by computing . 

In both cases the coefficients a, b, c are chosen such that 
the shock strength on the upper surface satisfies Equation (20) 
and these values are retained throughout the flow field. 


a 


(26) 
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5. NUMERICAL ALGORITHM 

The basic algorithm in the Ballhaus-Goor jian (ref. 5) 
computer code LTRAN2 solves the equation 



2M 


oo U 


^t = 0 


(27) 


in conservation form, where F(^ ) = [1 - M^ - (Y + l)M*^<j) ] . 

The present, modified, form of the TSD equation is in the 
form of Equation (27) [but with a different F ] and the same 
algorithm will suffice. However, because of the type dependent 
switching of the difference scheme it is important to ensure the 
conservation properties of the modified TSD solution. 


The general form of the type dependent differencing used in 
LTRAN2, which is applied only to the first term of Equation (27), 
is 


F 2 Ax 

X 


-1 


(F 


- F 


1 + 


^ “ 2 


l ) (1 - 


+ e..^(F 


1 - 


^ 2 -J 


(28) 

where is a switching operator given by 


0 if F^ > 0 

e . = s 

^ 1 if F. < 0 (29) 

t. 1 


The suffix i denotes the i^^ point and Ax is the mesh spacing. 
The differencing of Equation (28) is equivalent to 


which is a conservation form of the equation. In the present 
method this form of differencing is used for the term and 
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the values of (j)^ in F(<J)jj) are computed using central differencing. 
Otherwise the basic LTRAN2 algorithm is unchanged. 


6. DISCUSSION OF RESULTS 

As a first attempt at correcting the conventional unsteady 
small disturbance equation the correction was computed using 
Method I in which only the steady state solution is corrected. 

The first example is for the flow around a NACA 64A410 airfoil 
at M =0.72 and a = 2° + 2° sin kt, where k is the reduced 
frequency and is equal to 0.2. The steady result is shown in 
figure 2 and it can be seen that the corrected result is a 
considerable improvement on the conventional TSD result. There 
is some additional discrepancy at the leading edge due to the 
inadequate treatment of the boundary conditions in the thin 
airfoil approximation. 

In figure 3 the unsteady pressure distribution calculated 
by Option I of the present method is shown for two stations on 
the oscillatory cycle. It can be seen that the present method 
gives too large a shock excursion and that the shock strength 
in the foremost position is much too weak. This is probably 
a consequence of the correction being only valid for the steady 
state solution. The conventional small disturbance solution 
procedure diverges rapidly in the unsteady mode due to excessive 
shock motion. 

In figure 4 the steady pressure distributions around a 
NACA 0012 airfoil at = 0.8 and a = 1.25° for the present method, 
and the Euler equation solution of Sells (ref. 8) is shown. It 
can be seen that the present results agree satisfactorily with the 
Euler equation solution as regards shock location but that the 
upper surface pressures are too high and the lower surface pressure 
is too low. However, the improvement over the conventional TSD 
result is substantial. The disagreement between potential 
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equation results and Euler equation results in the leading edge 
region is also apparent in figure 4 if the results of Holst 
(ref. 9) are compared to the Euler equation results. Since TSD 
methods are usually "tuned" to approximate the full potential 
results it is possible that the discrepancy between the Euler 
solution and the present solution is due to fundamental differences 
in the numerical procedures used to solve full potential and 
Euler equations, since rotational effects should not be very 
important in the leading edge region. 

In figure 5 the unsteady pressure distribution around a 
NACA 64A410 airfoil at M- = 0.72 and a = 2° ± 2° sin kt at 

CO 

k = 0.2 is shown where Method II is used. It can be seen that 
this option considerably improves the shock behavior over that 
shown in figure 3. It should be noted that part of the 
difference in results of the present method and the Euler equation 
method is due to the different mesh size and in the computations 
this effects the shock capture properties of the algorithms. 

Finally in figure 6 a weak shock example is shown to 
illustrate the fact that the present method gives results close 
to the conventional TSD solution for weak shock waves. The 
airfoil is a NACA 0012 airfoil at M =0.8 and a = ± 1/2° sin kt 
with k = 0.2. 

Although the present method does give substantially improved 
results over the conventional TSD solution, there are still 
several points regarding both the shock details in unsteady 
motion and the pressure distribution in regions of the flow 
outside the shock regions. It is desirable to obtain further 
Euler equation solutions for such an investigation. 
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7. CONCLUDING REMARKS 


A theory to correct the Transonic Small Disturbance 
equation to treat strong shock waves in unsteady flow has been 
developed. Although a fairly complete thoery has been developed, 
only a simplified form has been computationally implemented. 

The comparisons of results of the present method with solution 
of the Euler equations is adequate as regards the shock location 
but in certain cases the pressure distribution elsewhere on the 
airfoil surface is not satisfactory. It is suggested that the 
discrepancy may be due to inherent differences in the numerical 
scheme used to solve both sets of equations. Finally, although 
the present method gives a considerable improvement over the 
conventional TSD theory, it is desirable to further test the 
present theory in order to fully evaluate the technique. 
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